##### SCRIPT QUE GERA A FIGURA 2 #####

# 1. CARREGA PACOTES ------------------------------------------------------

library(tidyverse)
library(janitor)
library(naniar)
library(mice)

# 2. CARREGA DADOS -------------------------------------------------------

#' Nesta etapa, foi implementada a leitura dos dados combinados de Enade-RAIS

# 3. IMPLEMENTA MODELOS DE REGRESSAO LOGISTICA ---------------------------

# 3.1 Modelo separado por área ----

imp_med <- filter(dados_graduados_empregados, area_curso_agreg == "Medicina")
imp_soc <- filter(dados_graduados_empregados, area_curso_agreg == "Ciências Sociais Aplicadas")
imp_ctem <- filter(dados_graduados_empregados, area_curso_agreg == "CTEM")

imp_dir <- filter(dados_graduados_empregados, area_curso_agreg == "Direito")
imp_edu <- filter(dados_graduados_empregados, area_curso_agreg == "Educação")
imp_eng <- filter(dados_graduados_empregados, area_curso_agreg == "Engenharia")

imp_hum <- filter(dados_graduados_empregados, area_curso_agreg == "Humanidades")
imp_sau <- filter(dados_graduados_empregados, area_curso_agreg == "Saúde e Bem-Estar")
imp_tec <- filter(dados_graduados_empregados, area_curso_agreg == "Tecnólogos")

f <- "ind_empregado ~ tp_sexo + raca_cor + idade_faixa + edu_fam + nt_ce_quartil_area + nt_fg_quartil + trab_grad + igc_faixa_4c + setor_ies_bin + regiao_curso"

model_med <- with(imp_med, glm(as.formula(f), family = binomial))
model_soc <- with(imp_soc, glm(as.formula(f), family = binomial))
model_ctem <- with(imp_ctem, glm(as.formula(f), family = binomial))
model_dir <- with(imp_dir, glm(as.formula(f), family = binomial))
model_edu <- with(imp_edu, glm(as.formula(f), family = binomial))
model_eng <- with(imp_eng, glm(as.formula(f), family = binomial))
model_hum <- with(imp_hum, glm(as.formula(f), family = binomial))
model_sau <- with(imp_sau, glm(as.formula(f), family = binomial))
model_tec <- with(imp_tec, glm(as.formula(f), family = binomial))

combine_med <- pool(model_med)
combine_soc <- pool(model_soc)
combine_ctem <- pool(model_ctem)
combine_dir <- pool(model_dir)
combine_edu <- pool(model_edu)
combine_eng <- pool(model_eng)
combine_hum <- pool(model_hum)
combine_sau <- pool(model_sau)
combine_tec <- pool(model_tec)

# 4. CALCULA PROBABILIDADES PREDITAS -------------------------------------

# 4.1 Por CE, area e setor da IES -----

# Funcao que cria dados-base para calculo de probabilidades preditas

prepara_new_data <- function(x = data, y) {
      x %>% 
            dplyr::distinct(area_curso_agreg) %>%
            filter(area_curso_agreg == y) %>% 
            rbind(dplyr::distinct(x, area_curso_agreg) %>% 
                        filter(area_curso_agreg == y)) %>%
            rbind(dplyr::distinct(x, area_curso_agreg) %>% 
                        filter(area_curso_agreg == y)) %>%
            rbind(dplyr::distinct(x, area_curso_agreg) %>% 
                        filter(area_curso_agreg == y)) %>%
            dplyr::arrange(area_curso_agreg) %>% 
            dplyr::mutate(tp_sexo = "Feminino",
                          raca_cor = "Negro/Indígena",
                          edu_fam = "Ensino Médio",
                          nt_ce_quartil_area = rep(c("q1", "q4"), 2),
                          nt_fg_quartil = "q3",
                          trab_grad = c("Não trabalha"),
                          igc_faixa_4c = "3",
                          idade_faixa = "25-",
                          setor_ies_bin = rep(c("Pública", "Privada"), each = 2),
                          regiao_curso = "SE") %>% 
            mutate_if(is.character, as.factor)
}

new_data_soc <- prepara_new_data(dados_graduados_empregados, "Ciências Sociais Aplicadas")
new_data_cte <- prepara_new_data(dados_graduados_empregados, "CTEM")
new_data_dir <- prepara_new_data(dados_graduados_empregados, "Direito")
new_data_edu <- prepara_new_data(dados_graduados_empregados, "Educação")
new_data_eng <- prepara_new_data(dados_graduados_empregados, "Engenharia")
new_data_hum <- prepara_new_data(dados_graduados_empregados, "Humanidades")
new_data_med <- prepara_new_data(dados_graduados_empregados, "Medicina")
new_data_sau <- prepara_new_data(dados_graduados_empregados, "Saúde e Bem-Estar")
new_data_tec <- prepara_new_data(dados_graduados_empregados, "Tecnólogos")

# Funcao que calcula probabilidades preditas

predict_prob <- function(x, y, z) {
      
      p <- x$analyses[[1]]
      p$coefficients = summary(y)$estimate
      predict(p, newdata = z,
              type = "response", se.fit = TRUE)
      
}

pred_soc <- predict_prob(model_soc, combine_soc, new_data_soc)
pred_cte <- predict_prob(model_ctem, combine_ctem, new_data_cte)
pred_dir <- predict_prob(model_dir, combine_dir, new_data_dir)
pred_edu <- predict_prob(model_edu, combine_edu, new_data_edu)
pred_eng <- predict_prob(model_eng, combine_eng, new_data_eng)
pred_hum <- predict_prob(model_hum, combine_hum, new_data_hum)
pred_med <- predict_prob(model_med, combine_med, new_data_med)
pred_sau <- predict_prob(model_sau, combine_sau, new_data_sau)
pred_tec <- predict_prob(model_tec, combine_tec, new_data_tec)

# Funcao que gera tabela com probabilidades preditas

prepara_tabela_2 <- function(x, y, z) {
      data.frame(z[,c("area_curso_agreg", "nt_ce_quartil_area", "setor_ies_bin")],
                 pred = x) %>%
            mutate(area_curso_agreg = y) %>% 
            select(1, 2, 3, 4, 5)
}

tab_trab <- bind_rows(
      prepara_tabela_2(pred_soc, "SOC", new_data_soc),
      prepara_tabela_2(pred_cte, "CTEM", new_data_cte),
      prepara_tabela_2(pred_dir, "DIR", new_data_dir),
      prepara_tabela_2(pred_edu, "EDU", new_data_edu),
      prepara_tabela_2(pred_eng, "ENG", new_data_eng),
      prepara_tabela_2(pred_hum, "HUM", new_data_hum),
      prepara_tabela_2(pred_med, "MED", new_data_med),
      prepara_tabela_2(pred_sau, "SAU", new_data_sau),
      prepara_tabela_2(pred_tec, "TEC", new_data_tec))

# 5. GERA GRAFICO --------------------------------------------------------

# Cria definicoes de layout do grafico

theme_plot <- theme(axis.text = element_text(colour = "black", size = 22),
                    axis.title.x = element_text(face = "bold"),
                    axis.title.y = element_text(face = "bold"),
                    legend.position = "bottom", 
                    legend.title = element_blank(), 
                    legend.text = element_text(size = 22),
                    text = element_text(size = 22),
                    strip.text = element_text(face = "bold", size = 22))

# Cria grafico

fig_2 <- tab_trab %>% 
      tibble() %>% 
      unite(setor_quartil, setor_ies_bin, nt_ce_quartil_area) %>%
      ggplot(aes(area_curso_agreg, pred.fit, fill = setor_quartil, color = setor_quartil)) +
      geom_point(stat = "identity", size = 4) +
      geom_errorbar(aes(ymin=pred.fit-pred.se.fit, ymax=pred.fit+pred.se.fit), width=.3) +
      scale_color_brewer(palette = "Set1") +
      scale_y_continuous(breaks = seq(from = 0.4, to = .9, by = .1), limits = c(0.4, .9)) +
      theme_bw() +
      coord_flip() +
      # scale_fill_grey(start = 0.8, end = 0.3) +
      labs(x = "", y = "Probabilidade predita de estar empregado formalmente") +
      theme(legend.position = "bottom", 
            legend.title = element_blank(),
            legend.text = element_text(size = 16),
            text = element_text(size = 18, color = "black"),
            axis.text = element_text(size = 18, color = "black"),
            axis.title = element_text(size = 18, color = "black")) +
      guides(colour = guide_legend(ncol = 2))

# 6. EXPORTA GRAFICO -----------------------------------------------------

jpeg(filename = "output/Figura 2.jpeg", width = 582, height = 388)
fig_2
dev.off()
